rm(list = ls())

# Parameters
W1.low <- 1
W1.high <- 1.2 
C.low <- 0.5
C.high <- 0.83
Overline.v <- 5

# Functions 
W0 <- function(v){
  out <- 1-pnorm(q = v, mean = 0, sd = 5)
  return(out) 
}

Admission <- function(v, 
                      overline.v = Overline.v,
                      mean = 2, sd = 2,
                      shape1 = 2, shape2 = 2){
  q <- (1-pnorm(overline.v, mean = mean, sd = sd)) / (1-pnorm(v, mean = mean, sd = sd))
  admission <- pbeta(q = q, shape1 = shape1, shape2 = shape2)
  return(admission)
}

Adj.stakes <- function(v,
                       w1 = W1.high, cost = C.high){
  adj.stakes <- cost / (w1 - W0(v))
  return(adj.stakes)
}

find.equilibrium <- function(w1 = W1.high, cost = C.high,
                             mean = 2, sd = 2,
                             overline.v = Overline.v,
                             shape1 = 2, shape2 = 2){
  f <- function(v){
    q <- (1-pnorm(overline.v, mean = mean, sd = sd)) / (1-pnorm(v, mean = mean, sd = sd))
    admission <- pbeta(q = q, shape1 = shape1, shape2 = shape2)
    adj.stakes <- cost / (w1 - W0(v))
    out <- admission - adj.stakes
    return(out)
  }
  lower <- 5*qnorm(1-w1+cost)
  upper <- overline.v
  out <- uniroot(f, lower = lower, upper = upper)$root
  return(out)
}

Pr.admitted <- function(apply, 
                        overline.v = Overline.v,
                        shape1 = 2, shape2 = 2,
                        mean = 2, sd = 2){
  admission <- Admission(v = apply, overline.v = overline.v,
                         mean = mean, sd = sd,
                         shape1 = shape1, shape2 = shape2)
  out <- (1-pnorm(q = apply, mean = mean, sd = sd))*admission
  return(out)
}

Pr.eligible <- function(apply, overline.v = Overline.v, 
                        mean = 2, sd = 2){
  out <- (1-pnorm(overline.v, mean = mean, sd = sd)) / (1-pnorm(q = apply, mean = mean, sd = sd))
  return(out)
}

find.critical.belief <- function(shape.1 = 2, shape.2 = 2){
  R <- function(p, 
                shape.p1 = shape.1, shape.p2 = shape.2){
    out <- pbeta(q = p, shape1 = shape.p1, shape2 = shape.p2) - 
      dbeta(x = p, shape1 = shape.p1, shape2 = shape.p2)*p
    return(out)
  }
  out <- uniroot(R, lower = 0.001, upper = 0.999)$root
  return(out)
}

Pr.apply <- function(v.hat, mean = 2, sd = 2){
  out <- 1-pnorm(v.hat, mean = mean, sd = sd)
  return(out)
}

V.2 <- seq(3.25, 5.5, length.out = 1000)


# Figure 
pdf('./figures/fig_2a.pdf', height = 7, width = 7)
plot(V.2, Pr.apply(v.hat = V.2), 
     type = 'l', ylim = c(0, 1), 
     lwd = 4, col = "#D4D4D4", xaxt = 'n', 
     ylab = 'Probability', 
     xlab = 'Threshold for applying', 
     cex.lab = 1.5, cex.main = 1.5)
axis(side = 1, at = c(find.equilibrium(w1 = W1.high), 
                      find.equilibrium(w1 = W1.low)), 
     labels = c(expression(paste(widehat(v)[P], '* ')), 
                expression(paste(widehat(v)[R], '* '))))
abline(v = find.equilibrium(w1 = W1.high), 
       col = 'black', lwd = 3, lty = 2)
abline(v = find.equilibrium(w1 = W1.low), 
       col = 'black', lwd = 3, lty = 2)
lines(V.2, Admission(v = V.2), 
      col = "#5E5E5E", lwd = 4)
lines(V.2, Pr.admitted(apply = V.2), 
      col = 'black', lwd = 4)
text('Permissive policy', 
      x = find.equilibrium(w1 = W1.high)-0.05, 
      y = 0.3, srt=-270)
text('Restrictive policy', 
      x = find.equilibrium(w1 = W1.low)-0.05, 
      y = 0.3, srt=-270)
legend(x = 3.2, y = 1, 
      col = c("#D4D4D4", "#5E5E5E", 'black'),
      legend = c('Pr(Apply)', 'Pr(Admission | Apply)', 
                 'Pr(Admission)'),
      bty = 'n', lwd = c(3, 3, 3))
dev.off()


pdf('./figures/fig_2b.pdf', height = 7, width = 7)
plot(V.2, Pr.apply(v.hat = V.2), 
     type = 'l', ylim = c(0, 1), 
     lwd = 4, col = "#D4D4D4", xaxt = 'n', 
     ylab = 'Probability', 
     xlab = 'Threshold for applying', 
     cex.lab = 1.5, cex.main = 1.5)
axis(side = 1, at = c(find.equilibrium(w1 = W1.high, cost = C.low), find.equilibrium(w1 = W1.low, cost = C.low)), 
     labels = c(expression(paste(hat(v)[P], '* ')), 
                expression(paste(hat(v)[R], '* '))))
abline(v = find.equilibrium(w1 = W1.high, cost = C.low), 
       lwd = 3, col = 'black', lty = 2)
abline(v = find.equilibrium(w1 = W1.low, cost = C.low), 
       lwd = 3, col = 'black', lty = 2)
lines(V.2, Admission(v = V.2), 
      col =  "#5E5E5E", lwd = 4)
lines(V.2, Pr.admitted(apply = V.2), 
      col = 'black', lwd = 4)
text('Permissive policy', 
      x = find.equilibrium(w1 = W1.high, cost = C.low)-0.05, 
      y = 0.3, srt=-270)
text('Restrictive policy', 
      x =  find.equilibrium(w1 = W1.low, cost = C.low)-0.05, 
      y = 0.3, srt=-270)
legend(x = 3.2, y = 1, col = c("#D4D4D4", "#5E5E5E", 'black'),
       legend = c('Pr(Apply)', 'Pr(Admission | Apply)', 'Pr(Admission)'),
       bty = 'n', lwd = c(3, 3, 3))
dev.off()
